Orientational ordering in hard rectangles: the role of three-body correlations 
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We investigate the effect of three-body correlations on the phase behavior of hard rectangle two- 
dimensional fluids. The third virial coefficient, B3, is incorporated via an equation of state that 
recovers scaled particle theory for parallel hard rectangles. This coefficient, a functional of the 
orientational distribution function, is calculated by Monte Carlo integration, using an accurate 
parameterized distribution function, for various particle aspect ratios in the range f — 25. A bi- 
furcation analysis of the free energy calculated from the obtained equation of state is applied to 
find the isotropic (I)-uniaxial nematic (N„) and isotropic-tetratic nematic (Nt) spinodals and to 
study the order of these phase transitions. We find that the relative stability of the N t phase with 
respect to the isotropic phase is enhanced by the introduction of B3. Finally, we have calculated the 
complete phase diagram using a variational procedure and compared the results with those obtained 
from scaled particle theory and with Monte Carlo simulations carried out for hard rectangles with 
various aspect ratios. The predictions of our proposed equation of state as regards the transition 
densities between the isotropic and orientationally ordered phases for small aspect ratios are in 
fair agreement with simulations. Also, the critical aspect ratio below which the N t phase becomes 
stable is predicted to increase due to three-body correlations, although the corresponding value is 
underestimated with respect to simulation. 

PACS numbers: 64.70.Md,64.75.+g,61.20.Gy 



I. INTRODUCTION 

The (two-dimensional) hard-rectangle (HR) model has 
recently received some attention due to the possibility 
that a dense film of these particles exhibits spontaneous 
tetratic order 0, 0, 0, 0- Additional interest originates 
from the fact that some types of organic molecular semi- 
conductors are made of rectangularly shaped molecules; 
a notable example is the PTCDA molecule, films of which 
have recently been studied quite intensely 5]. Even 
though the interactions between these molecules involve 
high-order polar forces (e.g. quadrupolar forces) it is of 
interest to investigate theoretically the intrinsic order as- 
sociated to purely excluded-volume effects with a view 
to predicting structural and thermodynamic properties 
of the film by incorporating other interactions via tradi- 
tional perturbation theories. The system we investigate 
in the present work mimics an incommensurate film of 
these molecules with only excluded-volume interactions 
involved and in the regime where molecules are free to 
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move in the film (i.e. fluid regime). Phases with two- 
dimensional crystalline order will be left for future work. 

Recently monolayers of various macroscopically-sized 
particles have been studied using mechanical vibrations 
on the monolayer to induce motion [|| . Even though this 
is an athermal, non-equilibrium system reaching steady- 
state configurations, these configurations are mainly 
driven by packing effects and should give the trend as 
to what types of order could be expected. In particular, 
tetratic order was observed in particles with sufficiently 
sharp corners, resembling rectangles, in contrast with 
particles such as discorectangles (projections of sphero- 
cylinders on the plane) which only exhibit nematic order- 
ing, or basmati-rice grains, which in addition may have 
a smectic phase. 

Apart from the possible interest in modelling the be- 
haviour of monolayers made of molecules with technolog- 
ical interest, our investigations have the additional, more 
fundamental aim of elucidating the effect that three- 
body correlations have on the orientational properties of 
two-dimensional fluids. Onsager showed that for three- 
dimensional hard rod fluids in the limit of infinite aspect 
ratio (hard-needle limit), n — > 00 (with n = L/a, L and 
a being the length and width of the constituent parti- 
cles), the ratio between the third virial coefficient and the 
second virial coefficient squared asymptotically vanishes, 
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for the isotropic fluid, as B 3 /B% ~ (a / L)\og(L / a) 0. 
Taking this result into account, he used a second-order 
virial expansion for the free energy as a functional of the 
orientational distribution function and obtained predic- 
tions for the isotropic (I)-nematic (N) phase-transition 
densities, exact in the above limit. By contrast, in two 
dimensions the above ratio between virial coefficients has 
the approximate limiting value of 0.514 [8J, implying that 
three-body correlations might play a very important role 
in the isotropic fluid even in the hard-needle limit; this 
is in sharp contrast with the three-dimensional case. An 
investigation of the effect of these high-order correlations 
on orientational ordering seems therefore appropriate. 

Since Onsager theory does not account for higher- 
than-two body correlations, an alternative theoretical ap- 
proach is required. Scaled particle theory (SPT), first 
developed for a mixture of hard spheres J3 and later ex- 
tended to anisotropic particles 0, 0, , includes as 
a main ingredient the exact analytic expression for the 
second virial coefficient 0, 0] , but again the third is 
approximated assuming Bj,fB\ — > in the hard- needle 
limit, an assumption that is incorrect. Also, it has been 
shown that, for a variety of particle shapes in two di- 
mensions, the fourth and fifth virial coefficients tend to 
negative values in the same limit 0, 0| . Thus it may 
very well occur that in these cases the virial series ex- 
hibits poor convergence, and a natural question arises: 
how does the phase behavior of anisotropic hard convex 
two-dimensional fluids change when three- and higher- 
body correlations, neither of which are included in the 
standard Onsager and SPT approaches, are taken into 
account? One of the aims of the present article is to shed 
some light on this question. For this purpose we develop 
an equation of state (EOS) for HR which exactly includes 
two- and three-body correlations in the nematic fluid and 
recovers SPT in the case of perfectly aligned particles. 

Recent investigations of the HR system have used the 
SPT approach 0,0 and Monte Carlo (MC) simulations 
0, to study its phase behavior. Aside from the usual 
isotropic- uniaxial nematic (N„) transition, these works 
have shown that this peculiar system exhibits a continu- 
ous transition between the isotropic phase and a tetratic 
nematic (N t ) phase. The latter is an orientationally or- 
dered phase but with D^h symmetry, i.e., the system is 
invariant under rotation of 7r/2. The SPT predicts that 
this phase is stable up to an aspect ratio k « 2.21 and 
that the packing fraction values of the I-Nj transition are 
around 0.85. This is in disagreement @ with the I-N t 
transition densities obtained from simulation for k = 1 
and 2 [3,|4j, which predicts values around 0.7. Consider- 
ing the importance that high-body correlations may have 
on the phase behavior of two-dimensional hard-convex 
bodies, it is the second purpose of this article to apply 
our model (which includes three-body correlations) to 
calculate the phase diagram of HR and compare the re- 
sults with those of SPT and MC simulations. The main 
features of the phase diagram are calculated using bifur- 
cation theory for the TN transition and also minimizing 



the nematic free energy functional resulting from our pro- 
posed EOS. 

The article is organized as follows. In Section II 
we describe the theoretical model for a general two- 
dimensional hard-convex fluid as applied to the isotropic 
and nematic fluids. Section III is devoted to the results 
obtained from the analysis of the theory, and compari- 
son is made with simulation results; also, the complete 
liquid-crystal phase diagram is presented. Finally, some 
conclusions are drawn in Section IV. The Appendix con- 
tains a detailed account of the bifurcation analysis and 
the minimization method used to analyse the phase be- 
haviour of the model, together with some details on the 
computer simulations. 



II. THEORY 

In the present section we introduce the theoretical for- 
malism necessary for the study of the phase behavior of 
the HR fluid. This includes the derivation of the EOS 
for the isotropic and nematic fluids, along with the cor- 
responding free energy density. The phase behaviour of 
the model, to be presented in the next section, is anal- 
ysed by means of two complementary techniques: a bi- 
furcation analysis of the free energy density, and a full 
minimization using an accurate functional form for the 
orientational distribution function. Details on these tech- 
niques are given in the Appendix. 



A. EOS for the isotropic fluid 

In this section an equation of state for the isotropic 
phase, to be extended later to the nematic phase, is 
proposed, on a somewhat ad- hoc, but at the same well- 
founded, basis. The virial coefficients are intimately re- 
lated to geometric properties of the planar objects mak- 
ing up the fluid. The second virial coefficient of planar 
hard particle has the analytic form 0] 

B 2 =v+^, (1) 

where v and C are the area and perimeter of the particle. 
Some approximate analytic expressions for third virial 
coefficients of isotropic fluids made of three-dimensional 
bodies, as a function of their volume, surface area and 
mean curvature, have been proposed |l5j |. When com- 
pared with results from numerical calculations [15j , some 
of these expressions are seen to constitute accurate ap- 
proximations. In two dimensions volume has to be sub- 
stituted by area, area by perimeter and mean curvature 
by a function proportional to the perimeter (the latter is 
true for the most representative two-dimensional convex 
bodies, i.e. rectangles, discorectangles, and ellipses). Fol- 
lowing some of the most successful approximations [l5| 
but translated to the two-dimensional case, we write the 
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following analytic expression for the third virial coeffi- 
cient 



B 3 = v 2 



47T 



-vC 2 



(4tt) 



-C 4 



(2) 



where the numerical coefficients a and (3 are chosen in 
such a way as to guarantee i) the correct asymptotic hard- 
needle limit, and ii) a good comparison with well-known 
EOS for some isotropic fluids, such as hard disks. For the 
latter we have C 2 — 4ttv so that the three terms in the 
right-hand side of Eqn. can be unified into the single 
term (1 + a + f3)v 2 . The SPT for hard disks is recovered 
by choosing a + (3 = 2, whereas the SPT form for B$ for 
a general anisotropic particle is obtained from J5J with 
a = 2 and ,8 = 0. Note, from Eqns. Q and J2J), that 
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in the infinite aspect-ratio limit, as the particle area 
is proportional to the product of the two characteris- 
tic lengths of the particle (the width a and the length 
L, the latter being the larger one), while the perimeter 
is proportional to their sum. Therefore we obtain the 
asymptotic limit B^/B 2 — > (3 when L/a — > 00 and a sen- 
sible choice is P — 0.514, the exact asymptotic value of 
this ratio |8|. 

The EOS is obtained by imposing two requirements: 
(i) the divergence of pressure at high packing fractions is 
of the form ~ (1 — r\)~ 2 as stated by SPT, and (ii) the 
second and third virial coefficients are obtained from the 
exact virial expansion 



0P = p + P 2 B 2 +p 3 B. 
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(where p is the density of particles). In other words, 
we require that the third-order virial expansion of the 
interaction part of the EOS, 



a 2 r] 2 + a 3 rj 3 

(i-v) 2 ' 



(5) 



(j] = pv being the packing fraction) coincides with the 
exact one @. This allows us to obtain a k (k = 1,2) 
as a2 = I + 62, and a 3 = b 3 — 2a 2 — 1, where the new 
coefficients 



B k 

b k = ^-l, k = 2,3, 
v 



(6) 



have been defined in terms of the virial coefficients B k . 
The resulting EOS has the form 



f3Pv 
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From this EOS the free energy density can be obtained: 
we first write 



(8) 



where tp = (fid+fex is the free energy per particle and c/?id 
and ip ex the corresponding ideal and excess contributions. 
Now using (3P from 0, Eqn. JBJ can be integrated to 
give 



0( V ) 



ln(l - 77) + 
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+ ln(l-77) 



b 2 + (b 3 -2b 2 )9( V ), (9) 
(10) 



The first two terms of Q and © are SPT-like terms. 
Eq. J7J), with the exact second and third virial coefficients 
for the particular case of parallel hard rectagles, recovers 
the SPT result [this is easily obtained if we substitute 
the exact values B 2 = 2v (b 2 = 1) and B3 = 3v 2 (63 = 2) 
in Eq. 0]. 

Inserting B 2 and the approximation for B3 from 
and J2J, respectively, in Eqns. Q and 10, we obtain our 
proposed EOS and the excess part of the free energy per 
particle for the isotropic fluid as 



(3Pv = 



n 



if 



1 — 77 (1 — 77)2 

— ln( 1 — 77) 1 1V 



1-7, 



7 [1 + (a -2 + &y)rj\, (11) 

+ 1 {a-2 + (3 1 )9(7 1 ), 

(12) 



where the anisometric parameter 7 = C 2 /(Awv) was de- 
fined. Note that for a — 2, f3 = 0, this equation recovers 
the SPT expression for hard convex bodies. From (jl 1|) 
the following expression for the reduced virial coefficients 
is obtained: 



Bt = 



B„ 



1 + [1 + {a - 1)(ti - 2)] 7 + (3(n - 2) 7 2 
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(13) 



B. EOS for the nematic fluid 



The EOS for the nematic fluid is now obtained from 
Eqn. J7J) by substituting the virial coefficients of the 
isotropic fluid B n (n — 2, 3) by their functional versions 
B n [h] for the nematic fluid; here h((f>) is the orientational 
distribution function. The latter coefficients are obtained 
from the definitions of B 2 and B3, in terms of integrals 
over the Mayer function: 



B k [h] = 



[I / Wthih) 



/C(0i 



K.((t>i,4> 2 ) 



dr/(r, <j> 12 ), 



(15) 



K{th.,<l>2,<h) = - ] d * J rfr7(r,^ 2 )/(r',^ 3 ) 
x/(r-r',0 13 ), (16) 

where tf> a p = <p a — is the relative angle between axes 
of particles a and /3, and /(r, <j> a f}) the Mayer function. 
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The corresponding free-energy functional ip[h] = ip\d[h] + 
(p ex [h] is obtained from Eqns. 10): 

tp^lh] = -ln(l - n) + -^—b 2 [h] + (b 3 [h] - 2b 2 [h}) 6{n) 

I '/ 

(17) 

with the ideal part exactly calculated from 

/■2tt 

<f id [h] = In?/ - 1 + / d(j)h{<j)) In [2irh(cj))} . (18) 
Jo 

The remaining virial coefficients are approximated from 
Eq. (0 by 

B n [h] = v n ~ 3 {(n - 2)B 3 [h] - (n - 3)B 2 [%} . (19) 

The integral over spatial variables in the definition of 
-B2[^] is known analytically for most convex bodies. In 
particular, for HR we have 

= (£ 2 + O I sin 0i2 1 + 2La (1 + | cos0 12 |) , 

(20) 

which is the excluded area between particles with relative 
orientation 0i2. However, the required double angular 
average over h((f>) has to be estimated numerically (we 
used Gaussian quadrature). Also, in the case of B 3 [h], 
all integrals have to be calculated numerically (using MC 
integration). For this purpose we found it convenient to 
use a parameterized orientational distribution function 

h((j>) = C exp X r cos(2r0), ^ (21) 

in terms of the n parameters X T (r = 1, . . . , n). C is a 
normalization constant. In practice two variational pa- 
rameters (n = 2) were used. Details on how this calcu- 
lation was realized in practice are relegated to the Ap- 
pendix. 




FIG. 1: Typical particle configurations as obtained from MC 
simulations, (a) isotropic phase; (b) tetratic phase; and (c) 
crystalline phase with particles arranged in one of the possible 
configurations. See text for details on the simulations. 



III. RESULTS 

In this section we present the main results obtained 
from the inclusion of three-body correlations into the 
EOS for the isotropic and the nematic fluids, as pro- 
posed in Section II. The results from bifurcation analysis 
and from numerical minimization of the free-energy func- 
tional are presented in Section III B. In the latter case we 
compare the results from the present theory with those 
obtained from SPT and from simulations. But before 
presenting the results, we show in Fig. a series of par- 
ticle snapshots extracted from our simulation runs. De- 
tails on the simulations are given later on. The three con- 
figurations are representative of an isotropic phase [Fig. 
0a)], a tetratic phase [Fig. 0b)], and a crystalline phase 
[Fig. 0c)]. 



A. Isotropic fluid 

Let us first compare the different approximations for 
B 3 and assess their quality according to the degree of 
agreement with MC-integration results for isotropic flu- 
ids made of different hard-convex bodies. Our own MC 
calculations have been carried out for hard rectangles, 
while those of Ref. 01 were focussed on hard discorect- 
angles and hard ellipses. All the results are plotted in 
Fig. [21 along with two different analytic approximations 
for the reduced virial coefficient B 3 . One of them (solid 
line) is calculated from Eqn. (|13|) with (3 — 0.514 [which 
gives the correct value of B 3 in the Onsager limit [l7| -see 
Eqn. (3)], and setting a = 1.611, which gives the correct 
third virial coefficient for hard disks (B^ w 3.125). The 
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other (dotted) line is also calculated from Eqn. (|13fl . but 
choosing a — 2 and (3 — 1/8, which reduces to the pro- 
posal made by Boublik in Ref . 0] . As will be shown be- 
low, this proposal approximates the EOS for the isotropic 
phase of hard convex bodies reasonably well. Also plot- 
ted in Fig. Elwith dashed lines are the best fits calculated 
from 

B 3 = Shbv + — — vC + j—r^C (22) 

47T [A-kY 

with values for the coefficients <5hb, ohb, and /3hb de- 
pending on the particle geometry (HB = HR, HDR, HE), 
i.e. hard rectangles, hard discorectangles, and hard el- 
lipses. Note that Eqn. Il22[| is the same as Eqn. J5J, but 
with a new numerical coefficient <5hb as a prefactor of v 2 . 
From Fig. [21 we can see that the present approximation 
(a = 1.611, j3 = 0.514) describes the behavior of B3 as a 
function of the anisometric parameter 7 much better than 
Boublik's proposal which, in the Onsager limit, gives the 
(wrong) value B% = 1/8. Also, if one is to describe the 
correct behavior of B% for different particle geometries in 
the whole range of 7, it is necessary to take (5hb 7^ 1- 

Numerical values for the coefficients B\ and B\ have 
been calculated in Ref. [l7j. for three different particle 
geometries, using MC integration. The results for B\ are 
shown in Fig. [31 Also plotted are our analytic proposal 
(solid line) and that of Boublik (dotted line). For HR 
with small anisometry values our approximation is bet- 
ter than that of Boublik, while the opposite occurs for 
high anisometries. For hard ellipses, Boublik's approach 
describes reasonably well the behavior of B\ in the whole 
range of 7 (except for very long particles which have neg- 
ative values of B%). 

The EOS obtained from the above approximations can 
be checked against MC simulations of systems of HR par- 
ticles. In order to realize this comparison we have car- 
ried out constant-pressure MC simulations on systems of 
HR with different aspect ratios in the range of pressures 
where the isotropic fluid is the stable phase. The results 
for k — 3 and 9 are shown in Figs. 01 (a) and (b), re- 
spectively. Simulations were done on systems of ~ 10 3 
particles, equilibrated along typically ~ 10 6 MC steps, 
and averaging over ~4x 10 6 MC steps. The system was 
prepared in each case in a crystalline low-density config- 
uration with perfectly aligned particles at low pressure. 
This configuration rapidly turned into a disordered con- 
figuration, which was then equilibrated. After averaging, 
the system was subject to a higher pressure and then 
equilibrated, and the process was repeated increasing the 
pressure. In this way the EOS in the entire region of 
isotropic stability was obtained. For comparison we also 
show in Figs. 0] (a) and (b) the EOS corresponding to 
the SPT (dashed line), Boublik proposal (dotted line), 
our proposal (solid line) and the EOS [Eqn. (7J)] with 
the virial coefficient B3 calculated from MC integration. 
As can be seen the SPT and Boublik's proposal approx- 
imate better the simulation results. However, given that 
both theories make wrong predictions of the behavior of 
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FIG. 2: Reduced virial coefficient _BJ as a function of the ani- 
sometric parameter 7. Simulation results are shown for hard 
rectangles (squares), discorectangles (asterisk), and ellipses 
(triangles). The solid and dotted lines are the results from 
Eq. (EH with (a,/3) = (1.611,0.514), and (a,/3) = (2,1/8) 
respectively. Also are shown with dashed lines the best fits 
from Eq. J23. 
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FIG. 3: Reduced virial coefficient B\ as a function of the 
anisometric parameter 7. All the lines and symbols label the 
same as in Fig. [21 



the third and fourth virial coefficients of HR as a func- 
tion of the anisometric parameter, these results are to 
be taken with caution in the sense that they could be a 
mere coincidence. We can also see from the figures that 
our proposal overestimates the pressure. This kind of be- 
havior is typical of fluids composed of hard-core particles 
which exhibit poorly convergent virial series; this seems 
to be the case for HR particles since their fourth and fifth 
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FIG. 4: Results from MC simulations (filled circles) on a sys- 
tem of 10800 HRs with k = 3 (a), and k = 9 (b). Dot-dashed 
line: Eq. Q with Bz calculated from MC integration, dashed 
line: SPT, dotted line: Boublik proposal, and solid line: our 
proposal [The B3 from Eq. © with {a = 1.611, /3 = 0.514}]. 

virial coefficients become negative for high anisometries. 

B. Bifurcation to nematic fluid 

We implemented numerically the bifurcation-theory 
analysis, described in detail in the Appendix, to calcu- 
late the spinodal instabilities from the isotropic phase 
to the N„ and Nt phases, and elucidated the order of 
these transitions within the same formalism. These re- 
sults were checked against a full minimisation of the free- 
energy functional employing the methodology outlined in 
Section IIB, which in addition enabled the N t — N u spin- 
odals, which cannot be easily calculated using bifurcation 
theory, to be obtained. Also, in order to have essentially 



exact results for the phase behaviour of this system, we 
performed constant-pressure MC simulations on systems 
of ~ 10 3 — 10 4 HR particles and obtained the equations 
of state and orientational order (details on these simula- 
tions are included in Section D of the Appendix). All of 
these results are described in the following. 

The I-N u and I-N t spinodal lines t]*(k) (the packing 
fraction at bifurcation as a function of the aspect ratio 
k) were calculated by solving Eqn. (39) for y = 77/ ( 1 — 77) 
(or 77) for a discrete set of values of k (see Appendix) . All 
the coefficients 5( fcl ' fc2 < fc3 ) that enter this equation were 
calculated via MC integration; typically ~ 10 s MC steps 
were used to evaluate these coefficients. To elucidate 
the order of transitions, we solved Eqn. (46) to find (i) 
the value ki at which the free energy difference between 
N u or Nt and isotropic phases changes from negative to 
positive, which in turn reflects the change of sign of the 
coefficient B* (see Appendix), and (ii) the value K2 for 
which the inverse of the isothermal compressibility of the 
N u or Nt phases [(>% 1 ' i; ) ] at the bifurcation point be- 
comes zero. Again, all the coefficients 63 (the rescaled 
third virial coefficient) and ^ fcl > fe2 ' fe3 ) ; necessary to solve 
Eqn. (46), were evaluated using MC integration with the 
same number of steps as previously. The quantities B* 
and (re^ are shown as a function of k in Figs. 4 (a) 
and (b) for the I-N u transition, and in Figs. 5 (a) and 
(b) for the I-N t transition. 

As can be seen from the figures, (x^ v) first changes 
sign from positive to negative at K2 ~ 4.62, and then 
diverges at k 4.11, coinciding with ki, the zero of B* 
[see Fig. 0(a)]. The latter has a pole at n — 3.23, which 
is the intersection point between the I-N u and I-Nt spin- 
odals (for larger values of k the I-N u spinodal lies below 
the I-N t spinodal). This result can be understood from 
the definition of B* [see Eq. I|37|l ]. which diverges at 
D* = 0; this in turn coincides with the condition A* = 
if we change k to 2k in Eqn. I|3()ll to include tetratic sym- 
metry. Thus B* as a function of k should diverge at the 
point where the I-N u and I-Nt spinodals intersect. The 
values of B* and v) as a function of k, calculated 
this time at the I-N t spinodal, are shown in Fig. H3 (a) 
and (b). From the figure we can see that they are al- 
ways positive, in particular for values of k less than 3.23. 
Thus we can conclude that the I-N t transition is always 
of second order. We also note the oscillatory behavior 
of B* as a function of k [see Fig. EJa)]; this feature has 
been shown not to be a consequence of numerical errors 
inherent to our MC integration: MC estimates of the co- 
efficients involved in the definition of B* were obtained 
by increasing the number of MC steps from 10 6 to 10 9 , 
and the oscillating behavior remained. 



C. Phase diagram 

The resulting spinodal instabilities from the I to the 
orientationally ordered phases, as calculated from the bi- 
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FIG. 5: The coefficient B* (a), and the inverse of the com- FIG. 6: The coefficient B* (a), and the inverse of the com- 



pressibility factor [x^ L vy (b) at the I-N u bifurcation point 
as a function of the aspect ratio calculated for a discrete set 
of values (open circles). The filled circles indicate the value 
of k for which they become zero. Thus, k* = K2 ~ 4.62 is the 
true tricritical point. 



pressibility factor [x^ L vY (b) at the I-Nt bifurcation point 
as a function of the aspect ratio calculated for a discrete set 
of values (open circles). Both magnitudes are always greater 
than zero, so the I-Nt is always of second order. The arrow 
indicates the maximum aspect ratio of Nt phase stability. 



furcation analysis, are shown in Fig. [7] Also shown in the 
same figure is the complete phase diagram resulting from 
SPT, already calculated in Ref. 0- As can be seen from 
the figure, the inclusion of three-body correlations con- 
siderably lowers the transition densities between isotropic 
and the orientationally ordered phases. The new results 
compare fairly well with those from MC simulations in 
the region of low particle aspect ratio (our simulations for 
k = 3, and simulations for k = 1 [3j and k = 2 repre- 
sented in the figure by open squares. Another interesting 
point to remark is that the critical value of k below which 
the N t phase is stable increases from n — 2.62 in SPT to 
K = 3.23 in the new theory. Finally, the I-N„ tricritical 
point occurs at n* = max(Ki,K 2 ) ~ 4.62 [see Fig. El (a) 
and (b)], which is lower than the SPT result (k* = 5.44). 
Thus, we can conclude that three-body correlations have 



the effect of lowering the transition densities, increasing 
the stability of the N t phase and making the I-N u transi- 
tion weaker. The enhanced stability of the tetratic phase 
is in agreement with simulation results. 

It is also apparent from Fig. 0that the I-N„ transition 
predicted by the new theory for high values of k occurs at 
packing fractions below those predicted by SPT. In the 
Onsager limit, the reduced transition density p r = p*L 2 
for the I-N„ transition can be calculated within the third 
virial-coefHcient approximation [see Eqn. I|47|l of the Ap- 
pendix]. This reduced density depends on the coefficient 
r* defined in Eqn. I|48l) . which can be calculated by ex- 
trapolating the data for t(k) obtained from simulations. 
These data, shown in Fig. [S] for high values of k, are 
fitted very accurately by means of a straight line that 
intersects the vertical axis at r* w 0.314. Inserting this 
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0.36 



FIG. 7: Phase diagram of the HR fluid. Continuous and 
first-order transitions are indicated by dashed and solid lines, 
respectively. Dotted lines indicate extension of I-N* line into 
region where Nt is preempted by uniaxial nematic phase. SPT 
transition lines and B$ spinodals are indicated by correspond- 
ing labels. Circles: minimisation of free-energy functional in 
B3 theory, giving first-order (filled circles) or second-order 
(open circles) transitions. Open squares: simulation results 
for the isotropic-to-nematic transition. 



value in @7J, we obtain p® 3 — 3.15, which is less than 
the SPT result pf. PT = 4.71 and much less than the MC 
simulation value, which has been estimated to be be- 
tween 7 and 7.5 3]. This disagreement is probably due 
to the poorly convergent character of the virial series. 
As already pointed out, the fourth and fifth virial coeffi- 
cients are negative in this limit, so the proper inclusion 
of higher-order virial coefficients is necessary in order to 
obtain an accurate approximation for the I-N transition 
densities. 

Another interesting aspect of the phase diagram is the 
failure of the new theory to reproduce the transition from 
the isotropic to the nematic phase in the range of large 
aspect ratios explored by our simulations. Note that, as 
will be discussed later, the simulations cannot reach any 
definite conclusion as to the real nature (whether tctratic 
or uniaxial) of the nematic phase, especially for large as- 
pect ratio. The fact that the isotropic-to-nematic tran- 
sition line rj{n) obtained from simulations in the range 
k = 1 — 9 is a smoothly decreasing monotonic curve and 
that this line is quite close to the spinodal line for the I- 
Nt transition obtained from the new theory in the same 
range of aspect ratios, may be indicating that the stabil- 
ity of the uniaxial phase is largely overestimated by the 
new theory, but that tetratic ordering is relatively well 
reproduced. This is simply a hypothesis not based on 
any real evidence. 



0.34 




0.32 



0.02 



FIG. 8: The coefficient t(k) [see Eq. I48|l ] as a function of k 
for a discrete set of k's (open circles) calculated via Monte- 
Carlo integration. The straight line calculated from mean 
square approximation intersects the ordinate at the value in- 
dicated in the figure. 



D. Further results 

In order to appreciate more deeply the differences be- 
tween the SPT and the new theory, we now compare the 
EOS for the isotropic and orientationally ordered phases 
and the behaviour of the order parameters, (71,92, with 
packing fraction. The latter are defined by 



dej) cos (2i(p)h((p), 



i = l,2 



(23) 



with qi the uniaxial order parameter and (72 the tetratic 
order parameter. The comparison is done in Figs. |5Ta- 
c) for the case k — 3 and in Figs. EH ( a_c ) for k = 9. 
Also, the MC simulation results are shown. In the case 
of the EOS, both theories severely overestimate the pres- 
sure in the nematic regime when k — 3; however, the 
transition point, as mentioned previously, is much better 
reproduced when three-body correlations are included in 
the theory. For the longer particles the pressures are 
better reproduced, but the location and nature of the 
transition from the isotropic to the nematic phase are 
not correct; as already mentioned, if the uniaxial ne- 
matic phase is not taken into account, three-body corre- 
lations seem to be very important in promoting tetratic 
order in the isotropic phase. These correlations alone, 
when higher-order correlations are not considered, prob- 
ably overemphasise the relative stability of the uniaxial 
nematic phase with respect to the tetratic phase in the 
case of long particles, and cause a premature instability 
of the latter as particles become longer. 

Comparison of the orientational distribution functions 
in the case k — 3 indicates again the role of three-body 
correlations. Fig. 1101 shows the corresponding func- 
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tion for the uniaxial nematic phase that coexists with 
the tetratic (new theory) or isotropic (SPT) phase; even 
though the new theory predicts a much lower transition 
density than SPT, tetratic ordering is much more pro- 
nounced in the new theory since at this value of k the 
tetratic phase is still stable. 

A point worth mentioning is the identification of the 
value of aspect ratio where the tetratic phase is no longer 
stable. The nonequilibrium macroscopic experiments by 
Narayan et al. .6] find substantial tetratic correlations 
in cylindrical particles with aspect ratio k = 12.6. Our 
present simulation data are not sufficiently detailed to 
give conclusive results. However, data for k = 7 (not 
shown) and k = 9 seem to be compatible with N t sta- 
bility: the value of the uniaxial order parameter q\ is 
compatible with zero in the whole density range studied. 
In the case k — 9 [Fig. HUf b)]. however, there seems 
to be some tendency in the uniaxial order parameter to 
increase from zero. 

Nevertheless, it is very difficult within our present anal- 
ysis, to settle this question. It is in fact difficult to dis- 
tinguish between the N t and N„ nematic phases, since 
the latter exhibits substantial tetratic correlations even 
for the longer particles considered (k = 9). Before fur- 
ther work is undertaken, all we can say conclusively from 
the simulation results is that at some packing fraction 
the isotropic phase begins to display substantial tetratic 
order in a rather abrupt manner; whether this order cor- 
responds to uniaxial or strictly tetratic nematic phases is 
a matter that would require more simulation work using 
e.g. larger systems. Our limited study is only intended 
to provide approximate phase boundaries for systems of 
particles with various aspect ratios (note that previous 
simulation work on this and related systems 0, 0, 
were more detailed, but restricted to a particular aspect 
ratio). 



E. Dependence on EOS adopted 

It should be noted that the packing fraction values 
of the I-N„ it phase transitions depend sensitively on the 
approximation used for the EOS of the HR fluid. This 
can be easily shown if we approximate the third virial 
coefficient as a function of the second virial coefficient 
using the relations 
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l = b 2 , Pj 2 = b 3 -ab 2 , (24) 

which can be easily obtained from and (J2J and are 
only valid for the isotropic fluid. From (|24[l we obtain 

b 3 =ab 2 +0bl (25) 

which can be used as an approximation of the third virial 
coefficient of the nematic fluid. Thus, inserting the above 
expression in Eq. Q, and carrying out the bifurcation 



FIG. 9: Results from the SPT, the present model and com- 
puter simulation for HR fluid with n = 3. (a) equation of 
state resulting from SPT (dashed line), our proposal (solid 
line) and MC simulation results (filled circles) . The arrow in- 
dicates the packing fraction of the I-N transition as estimated 
by simulation, while the filled square indicates the location 
of the I-N bifurcation point predicted by the present model; 
(b) behaviour of the uniaxial (solid line) and tetratic (dashed 
line) order parameters with packing fraction. Results from 
SPT and our model are indicated by the corresponding label. 
Symbols are simulation results for the order parameters (open 
symbols: uniaxial, filled symbols: tetratic); (c) orientational 
distribution functions at the coexistence packing fraction for 
the uniaxial nematic phase, from SPT (dashed line) and the 
present model (solid line). 
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<|>/Jt 

FIG. 10: Same as in Fig. but for k = 9. In (c) are shown 
the orientational distribution functions calculated at packing 
fractions separated from the bifurcation packing fraction a 
relative distance (A = 0.01467) equal to that between the 
isotropic and nematic coexisting packing fractions for k = 3. 
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FIG. 11: Orientational distribution function h((f>) from sim- 
ulation for the case k — 9 and packing fraction rj = 0.618. 
Symbols: simulation data. Line: best fit to Eqn. I|49|l . Re- 
sulting values for the order parameters are qi = 0.153 and 
q 2 = 0.628. 



analysis described in Section IIC, we arrive at 



hi 



4 

*(») = V 



k = 1,2 



2(3 



(26) 
(27) 



for the second order expansion of the free-energy differ- 
ence between the I and N phases at the bifurcation point; 
here the subindex k = 1, 2 labels the N u and N t phases, 
respectively. Solving equation Aip^ = for ^(k), we 
find the spinodals shown in Fig. El for different values 
of (a, (3) corresponding to those of SPT, Boublik's pro- 
posal, and our proposal for the EOS of the isotropic fluid. 
In the same figure the spinodal line resulting from the 
EOS with the exact third virial coefficient is also plot- 
ted. Comparing the latter with those obtained from the 
different approximations embodied in (|25l) . we conclude 
that the location of the I-N u tricritical point changes only 
if one uses the exact three-body correlations. The reason 
for this behavior can be elucidated from Eqn. Ij26(l : the 
values of (77* , «;*) calculated from A<^i = At/32 = gives 
us n* = (3 + s/5)/2, independent on the choice of (a, (3) 
as the function ^(y) does not depends on k. 



IV. CONCLUSIONS 

The main results presented in this article can be sum- 
marized as follows, (i) The inclusion of many- (higher- 
than-two) body correlations in two-dimensional systems 
of hard anisotropic bodies is of crucial importance in or- 
der to adequately describe the phase behavior of these 
systems. In two-dimensions two-body interactions are 
not enough to make quantitative predictions of their 
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FIG. 12: Spinodals of the transitions between the isotropic 
and orientational ordered phases. From top to bottom are 
shown the results from the SPT, the EOS with the approxi- 
mated 63 (the Boublik, and our proposal), and from the EOS 
with the exact 63. The dotted lines represent the position of 
the tricritical points common to all the approximations and 
that corresponding to the exact B3. Also are shown the sim- 
ulation results. 



covers SPT in the limit of spatially uniform phases and 
combines Onsager and fundamental-measure theories, in- 
dicate that the tetratic phase is preempted (in the sense 
of bifurcation theory, i.e. spinodal lines) by a spatially 
ordered phase. Inclusion of three-body correlations could 
severely affect this result since these correlations may af- 
fect both phases differently. In fact, simulations avail- 
able so far support the conclusion that the tetratic phase 
may be stabilised prior to crystallisation. However, even 
in the case that it were possible to construct a density 
functional, suitable for such non-uniform phases, and in- 
corporating three-body correlations, the effort involved 
in the minimization with respect to the full density pro- 
file p(r, (j>) would be rather huge. Work along this avenue 
is now in progress in our group. 



Acknowledgments 

Y.M.-R. was supported by a Ramon y Cajal research 
contract from the Ministerio de Ciencia y Tecnologfa 
(Spain). This work is part of the research Projects No. 
BFM2003-0180, FIS2005-05243-C02-01, FIS2005-05243- 
C02-02 and FIS2004-05035-C03-02 of the Ministerio de 
Education y Ciencia (Spain), and S-0505/ESP-0299 of 
Comunidad Autonoma de Madrid. 



phase behavior, a crucial difference with respect to three- 
dimensional systems. This conclusion is supported by 
simulation results, ii) We have proposed an EOS and 
a corresponding free energy density functional for fluids 
of hard rectangles that incorporates three-body correla- 
tions. While predicting pressures for the isotropic and 
nematic fluids which are too high when compared with 
simulation values, the theory gives values for the coexis- 
tence densities of the I-N t transition that compare fairly 
well with the simulation results for small values of k. A 
shortcoming of the theory is that the third virial coeffi- 
cient, which is incorporated exactly, has to be evaluated 
numerically beforehand. This is a practical, not funda- 
mental, limitation of the theory, which can be circum- 
vented in all cases (i.e. for all different particle geome- 
tries in two dimensions). 

A striking prediction of the theory is the stability of a 
tetratic phase, in good quantitative agreement with sim- 
ulations for low particle aspect ratio. We conclude from 
this that the four-fold correlations present in this phase 
are basically taken care of by our third-virial coefficient 
based theory, but not by the usual scaled-particle theory, 
which incorporates only two-body correlations. More at 
variance with simulation results is the case of high aspect 
ratios. In this regime the stability of the uniaxial nematic 
phase is overestimated with respect to the isotropic fluid. 

As a final comment, we must say that no attention 
has been paid to non-uniform phases (smectic, columnar 
or solid) in the present work. Previous studies by our 
group |2j , based on a density- functional theory which re- 



V. APPENDIX 

A. Bifurcation analysis 

In this section we introduce the formalism that allowed 
us to calculate the spinodal instabilities and the order of 
the phase transitions between the isotropic and orienta- 
tional ordered phases. This formalism is quite general, in 
the sense that it is independent of the geometry of parti- 
cles, and includes two- and three-particle correlations. As 
usual, the analysis starts from an order-parameter expan- 
sion of the free-energy difference (A<p) between the bifur- 
cated (orientationally ordered) and the parent (isotropic) 
phases about the bifurcation point, and then the evalu- 
ation of the inverse isothermal compressibility (>f~ 1 ) of 
the bifurcated phase at the same point. The existence 
of a tricritical point, at which the order of the transi- 
tion changes from second to first order, is predicted from 
the first change of sign (from positive to negative) of Atp 
or k^ 1 . The starting point of the bifurcation analysis 
is to assume that the orientational distribution function 
near the I-N bifurcation point can be approximated as a 
Fourier series in small amplitudes hk ~ e fc (where e is an 
small parameter), truncated at second order, i.e. 

h{<p) ps -(l + /iicos20 + /i2cos40). (28) 

Inserting this expression into 10 and (|18|l , we obtain the 
difference between the nematic and isotropic free energies 
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per particle as 

Aip = ip N - ipi w Ah\ + Ch\h 2 + Dhj + Ehf,(29) 
where the coefficients A, C, D. E have the form 



l r 

= 4 
1 
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(2,2) 



(30) 

(31) 
(32) 
(33) 



and we have defined y = rj/(l — rf). 6(y) — y — ln(l + y) 
is the same function as in (|10fl . but here in terms of the 
new variable y. The coefficients 



Y[ / d(f>i cos 2ki4>i 



i), i = 2,3, /cjGN,(34) 



have also been defined, which originate from two- (i = 2) 
and three- (i = 3) body correlations. For i = 2 one can 
obtain analytic results which, for the specific case of HR, 



give 



, (fe.fe) _ 



(L + (-l) k a) 



(4k 2 - 1)tt 



(35) 



If we set 6{y) = in the SPT result is recovered. 

Minimizing the free energy difference l|29l) with respect to 
h 2 , we obtain h 2 as a function of hi [h 2 = —Ch\j (2D)] 
which, inserted in (|29|l . results in 



Aip 
B 



Ahj+Bhj, 
4D 



E 



(36) 
(37) 



where A and B are functions of the variable y. Minimiz- 
ing Eq. (|36[) with respect to hi, and taking into account 
the expansion of y about its bifurcation value y* , i.e. 
y w y* + y( 2 'hf, we arrive at 



2/ii 



A* + (2B* + A* y yW) hf \ = 0, (38) 



where A y is the first derivative of A with respect to y, and 
the asterisk on A, B, and A y means that these functions 
are evaluated at the bifurcation point y* . Solving l|38|) 
order by order, we obtain two equations: 



A* = 0, 



y 



,'2) _ 



2B 



(39) 
(40) 



the first one allowing to find y* , and hence the packing 
fraction T)*(k), as a function of the particle aspect ratio 



K = L/a, i.e the spinodal line of the I-N phase transition. 
Expanding (|36[) about the bifurcation point, and using 
||s2]J and gOJl, we obtain 



Aip = -B*h\, 



(41) 



which indicates that the I-N transition is of first order if 
B* < 0. Eqn. (|41|) can be written in a different, more 
convenient form, with use of h\ = (y — y*) /y^ and Eqn. 
(|4tjp . which results in 



(42) 



Using the definition of the inverse isothermal compress- 
ibility, x^ 1 = pd((3P)/dp, in terms of the y variable, 



■ v = y(l + y)— [ y — 



dy 



together with Eqn. I)42|) . we find 



(X^V)* = (xf 1 «)* -(y*) 3 (l + y*) 



2B* 



(43) 



(44) 



where 



(k^v)* = y*(l + y*)(l + 2y*b 2 ) 



(y*f(3 + 2y*) 



l 



y 



(b 3 - 26 2 ). (45) 



The existence of a tricritical point, at which the I-N 
transition changes from second to first order as particles 
change from large to small aspect ratios, can be found for 
a value of the aspect ratio k* satifying k* — max (k 1; k 2 ), 
where Kj (j — 1, 2) are the solutions to the equations 



S*(«i)=0, 



fas* 1 ") 



«2 







(46) 



The preceding analysis with hk ^ (k = 1,2) corre- 
sponds to the bifurcation analysis of the transition be- 
tween the isotropic and the uniaxial nematic phase N„. 
If hi = 0, h 2 ^ 0, the bifurcating phase is a tetratic ne- 
matic phase Nf. To carry out the bifurcation analysis 
for the I-N t transition, we can use exactly the same for- 
malism, except that we have to make the substitutions 



'2A-- 



and b 



(ki,...,ki) 



,(2ki 



,2ki 



in Eqns. lt2D|l - 



hk - 

Taking the Onsager limit k — > 00 in Eqn. (|30|) and 
considering that, in the asymptotic limit [see Eqn. ©], 

the coefficients b^'^ and 63 1 ' ' are of order k and n 2 , 
respectively, the condition l|39|) is equivalent to solving a 
second-order equation with respect to the reduced den- 
sity p r = p* L 2 , with the solution 



1 



p r 



(VI + 3ttt* - 1) , 



where we have defined the coefficient 



lim t(k), t(k) 



37T b 



(1,1,0) 



(«) 



(47) 



(48) 
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The limit t* — > of Eq. igJJ recovers the SPT result 
p r = 37r/2. Also p r (r*) as a function of t* is a mono- 
tonically decreasing function whose domain and image 
arc [— 1/37T, oo) and (0, 3tt], respectively. Thus if r* > 
(r* < 0), the I-N transition in a two-dimensional hard- 
needle fluid occurs at a reduced density in the interval 
(0,3tt/2] ([37r/2,37r)). 

B. Calculation of B3 

The third virial coefficient B 3 ({X T }) was obtained by 
MC integration using, for the orientational distribution 
function, the form 

h(4>) = Cexp (Ai cos 20 + A 2 cos 40) (49) 

which contains two free parameters, Ai and A2. The MC 
data for B3 were obtained for each value of k explored 
and for fixed values of Ai and A2. The technique fol- 
lowed was a generalization of the standard method for 
isotropic fluids jl9j : each step involved generating an- 
gles for the three rectangles [see Eqn. (II tip ] and positions 
for two of them [the first rectangle is placed at the ori- 
gin, see Eqn. (|16|) ]. Since the B 3 coefficient involves a 
single irreducible cluster integral where all three rectan- 
gles overlap, the positions of the second and third rectan- 
gles were chosen within their respective excluded volumes 
with the first rectangle to insure overlap, and only one 
overlap condition (second with third rectangles) had to 
be checked. Angles were generated using an acceptance- 
rejection method according to the angular distribution 
function h((j>) corresponding to the values of Ai and A2 
(the method was checked by computing the second virial 
coefficient B2 for the isotropic case, which can be com- 
pared with the exact result, and also the third virial co- 
efficient for some special particle orientations where this 
coefficient is analytic; in this case the computed values 
of B2 and B 3 for high values of Ai and A 2 = tended to 
the correct value). A single MC step involves generating 
one chain of rectangles, and 1 — 2 x 10 7 MC steps were 
used in the calculations for each set of values of Ai, A 2 . 

We generated numerical values of B 3 at a collection 
of mesh points on a rectangular region of the Ai — A2 
plane. The extension of this region was chosen according 
to the values of the packing fraction. In any case it con- 
tained the origin (Ai = A2 = 0) to allow for the isotropic 
phase. In some cases the region [—1,1] x [—1,1] did suf- 
fice; in others, higher minimum and maximum values for 
the parameters were needed, especially when a first-order 
transition was detected. The mesh interval was typically 
A A = 0.1, with finer meshes when required. In order 
to use these data in a practical way, the data were fit- 
ted in two different ways. One involves constructing a 
polynomial P N (Xi, A 2 ) = £n=o ELo c nm Xf\^~ m by a 
least-square procedure. Symmetry considerations require 
some terms of this polynomial expansion not to appear, 
and the corresponding coefficients c nm were taken to be 
zero. The degree of the polynomial was typically in the 



range 8 — 10. In the case of the I-N t transition, which 
only involves 92 (and hence A2), calculations were also 
done using fits to a polynomial depending only on A2 
(since necessarily Ai = 0). Results are consistent with 
the previous results based on a full fitting. 

For high packing fractions a fit to a function in the or- 
der parameters (91, 92) is more suitable since their values 
are close to one, whereas the A parameters grow with- 
out limit. However, the dependence of B 3 on ((71,92) is 
strong. We found it useful to use a combination of poly- 
nomials in the 9's and factors of the form (<fc ± 1)™, with 
n a power whose value is optimized in the fit. 



C. Minimization of the free-energy functional 

The minimizations were done using a variational 
scheme. An important question is how accurate is the 
variational function Q49J1. We can assess the quality of 
this function by comparing with simulation results for the 
distribution function h((f>). Fig. II II shows a distribution- 
function histogram obtained from a constant-pressure 
MC simulation, over 2 x 10 6 steps, of a fluid with k = 5 
at pressure Pa 2 /kT = 1.4. This is clearly a nematic 
phase with tetratic order (whether this corresponds to 
a uniaxial or purely tetratic phase is a different mat- 
ter; extremely long runs are probably needed to fully 
equilibrate the system. For the present purpose this is 
of no importance). A least-square fitting to the varia- 
tional function gives Ai = 0.033, A 2 = 1.217 (q x = 0.008, 
92 = 0.523), which results in the function represented in 
the figure. Histograms exhibiting more structured ori- 
entational order can be similarly fitted with comparable 
accuracy. The function 14911 is therefore suitable as a 
variational function. 

The nematic order parameters can be related to the 
variational parameters via Eqn. (|23l) . There is a 
one-to-one correspondence between the sets (91,92) and 
(Ai,A 2 ). In our calculations the function ip(qi,q2;K,rj) 
was then minimized with respect to 91 , 92 , using a stan- 
dard Newton-Raphson technique. All transitions are ob- 
tained as second-order transitions, except in the approx- 
imate interval 3 < k < 5 where discontinuous transitions 
were found. Of course these results are consistent with 
those from bifurcation theory, which is otherwise better 
suited for the calculation of the tricritical points since 
it is not tied to any variational scheme. The value of 
the functional minimization can be better appreciated in 
the case of the N t -N u transition, which cannot easily be 
obtained using bifurcation theory. 

The results of the minimisation using the resulting fit- 
ted function for B 3 are very sensitive to details such as 
mesh interval and degree of fitting polynomial. A de- 
tailed study of the whole procedure, including fine tuning 
of the above parameters, was therefore necessary. The 
accuracy of the results is sufficient to locate the phase 
transitions with respect to packing fraction, and to dis- 
criminate between first- and second-order phase transi- 
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tions when the system is far from the tricritical points, 
but the exact density gap in first-order transitions (which 
is otherwise small) could not be obtained with the present 
numerical implementation. 

D. Some details on MC simulations 

Our constant-pressure MC simulations were performed 
on systems of ~ 10 3 — 10 4 HR particles, using rectangu- 
lar cells and periodic boundary conditions. The equa- 
tion of state, orientational distribution function and ne- 
matic order parameters were obtained during the course 
of these simulations. The simulations were run typically 
over ~ 10 7 MC steps for equilibration and ~ 2 x 10 7 
MC steps for averaging (slow orientational dynamics, es- 
pecially in the case of long particles, require longer runs 
than in the isotropic phase). The value of the pressure 
is fixed at some constant value. The average density (or 
packing fraction 77) is obtained, which gives the equation 
of state P(r]). 

The orientational order is obtained from the eigenval- 
ues and eigenvectors of the order matrix, defined for a 
given particle configuration as 

1 N 

s « = ]vE( 2 ".- #i *-*«) ( 5 °) 
fe=i 

where n fc is a unit vector along the long axis of the fcth 
particle. These is to be averaged over MC configura- 
tions. The eigenvector associated with the largest eigen- 
value gives the direction of the primary director, and this 
eigenvalue is the uniaxial order parameter, q±, which can 
also be calculated from the average 

1 N 

k=l 

where <po is the polar angle of the director (which depends 
on the configuration) and <f>k the polar angle of the parti- 
cle long axis, both with respect to some fixed direction in 
the plane. The tetratic order parameter can now be cal- 
culated from a similar equation, with the factor 2 in the 



cosine substituted by 4. Also, the orientational distribu- 
tion function h(<fi) was calculated as a histogram, using 
the angle 4>q as the origin. From this the order parame- 
ters qi can also be obtained from Eqns. (|23|) . Yet another 
route is provided by the asymptotic value of the orien- 
tational correlation functions. No attempt was made at 
calculating these functions in the present work. 

The function h(<fi) is always seen to have two max- 
ima: a primary and a secondary maximum, separated by 
7r/2; in the uniaxial nematic phase these maxima should 
have different heights, whereas in the tetratic phase their 
heights should be statistically equal. Due to many ef- 
fects that affect the simulations, distinguishing these two 
situations is a rather delicate problem and our limited 
study did not allow identification of the true nature of 
the nematic phase. Questions such as effect of boundary 
conditions, system size, etc, may be of paramount im- 
portance in this analysis. For example, the rectangular 
periodic boundary conditions used in this work could be 
artificially promoting tetratic ordering in the system. 

As is characteristic of two-dimensional systems with 
continuous symmetries, the orientationally ordered 
phases of the present model seem to exhibit quasi-long- 
range order at long distances 0. This means, in par- 
ticular, that the order parameters may show a strong 
system-size dependence. This point has not been consid- 
ered at all, since our only aim was to establish approxi- 
mate phase-stability boundaries that could serve as a test 
bed against which the (otherwise approximate) theories 
could be tested. 

Finally, since this was not the aim of this work, and 
also due to the difficulties involved in dealing with pos- 
sibly multiply degenerate structures, both periodic and 
nonperiodic, crystalline configurations have not been 
studied in any detail; however, freezing into glassy states 
in compression runs were observed to occur (these states 
were characterised by extremely low particle diffusion) 
at high density. These densities at quite close to those 
at which melting into a nematic phase from crystalline 
configurations are observed to occur in expansion runs 
(starting from crystals with various types of packing - 
particles perfectly aligned on a rectangular lattice, square 
clusters on square lattices, various random tilings, etc.) 
A full discussion of this issue can be found in Ref . [4( . 
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